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We present a new class of exponential integrators for ordinary dif- 
ferential equations: locally exact modifications of known numerical 
schemes. Local exactness means that they preserve the linearization 
of the original system at every point. In particular, locally exact in- 
l/^ . tegrators preserve all fixed points and are A-stable. We apply this 

On I approach to popular schemes including Euler schemes, implicit mid- 

-vj ' point rule and trapezoidal rule. We found locally exact modifications 

■' ■ of discrete gradient schemes (for symmetric discrete gradients and co- 

f^ , ordinate increment discrete gradients) preserving their main geomet- 

ry ' ric property: exact conservation of the energy integral (for arbitrary 

multidimensional Hamiltonian systems in canonical coordinates). Nu- 
merical experiments for a 2-dimensional anharmonic oscillator show 
^ I that locally exact schemes have very good accuracy in the neighbour- 

H ' hood of stable equilibrium, much higher than suggested by the order 

■ ■ ■ of new schemes (locally exact modification sometimes increases the 

order but in many cases leaves it unchanged). 
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1 Introduction 

The motivation for introducing "locally exact" discretizations is quite natu- 
ral. Considering a numerical scheme for a dynamical system with periodic 
solutions (for instance: a nonlinear pendulum) we may ask whether the nu- 
merical scheme recovers the period of small oscillations. For a fixed finite 
time step the answer is usualy negative. However, many numerical scheme 
(perhaps all of them) admit modifications which preserve the period of small 
oscillations for a fixed (not necessarily small) time step h. An unusual fea- 
ture of our approach is that instead of taking the limit /i — )■ 0, we consider 
the limit x^ ~ a; (where x is the stable equilibrium). As the next step we 
consider a linearization around any fixed x (then, in order to preserve the 
condition Xn ~ 2;, some evolution of x is necessary). In other words, we com- 
bine two known procedures: the approximation of nonlinear systems by linear 
equations and explicit exact discretizations of linear equations with constant 
coefficients. An essential novelty consists in applying these procedures to a 
modified numerical scheme containing free functional parameters. The case 
of small oscillations was presented in [8j. Our method works perfectly for 
discrete gradient schemes [201 1221 |29]. We succeded in modifying discrete 
gradient schemes in a locally exact way without spoiling their main geomet- 
ric property: the exact conservation of the energy integral. The preservation 
of geometric properties by numerical algorithms is of considerable advantage 
[HllTT]. Promising results on one- dimensional Hamiltonian systems can be 
found in [9l [10] (by one- dimensional Hamiltonian system we mean a Hamil- 
tonian system with one degree of freedom). In this paper we extend our 
approach on the case of multidimensional canonical Hamiltonian systems. 
Moreover, we present locally exact modifications of forward and backward 
Euler schemes, implicit midpoint rule and trapezoidal rule. 

A notion closely related to our local exactness has been proposed a 
long time ago [27], see also [21]. Recently, similar concept appeared under 
the name of linearization- preserving preprocessing [23] , see also below (sec- 
tioning])- Our approach has also some similarities with the Mickens approach 
[24j, Gautschi-type methods [1^ and, most of all, with the exponential inte- 
grators technique [3l[T5l[25]. The definition of an exponential integrator is so 
wide (e.g., "a numerical method which involves an exponential function of the 
Jacobian", [15]) that our schemes can be considered as special exponential 
integrators. In spite of some similarities and overlappings, our approach dif- 
fers from all methods mentioned above. In particular, according to our best 



knowledge, discrete gradient schemes have never been treated or modified 
in the framework of exponential integrators and/or linearization-preserving 
preprocessing. 

2 Locally exact numerical schemes 

We consider an ordinary differential equation (ODE) with the general so- 
lution xit) (satisfying the initial condition x(to) = Xq), and a difference 
equation with the general solution Xn- The difference equation is the exact 
discretization of the considered ODE if Xn = x{tn). 

2.1 Exact discretization of linear systems 

It is well known that any linear ODE with constant coefficients admits the 
exact discretization in an explicit form [28] , see also [H [TJ [23] . We summarize 
these results as follows, compare [10] (Theorem 3.1). 

Proposition 2.1. Any linear equation with constant coefficients, represented 
in the matrix form by 

I ^^.6. (2.1) 

(where x = x{t) G M'^, b = const e M.'^ and A is a constant invertible real 
d X d matrix) admits the exact discretization given by 

Xn+i - Xn = {e^-^ - I)A-^ {AXn + 6) , (2.2) 

where hn = tn+i — tn is the time step and I is the identity matrix. 

Corollary 2.2. We may look at the exact discretization Ii2.3\) as a modifica- 
tion of the forward Euler scheme. Indeed, we can rewrite ^2.2) as 

K,^ (x„+i - Xn) = AXn + b , (2.3) 

where A„, defined by 

A„ = A-i(e'^"^-/) , (2.4) 

is a matrix "perturbation" of the time step hn. Indeed, A„ = /i„/ + 0{h^. 



Exact discretizations found direct application in numerical treatment of 
the classical Kepler problem |31 El E]- The exact discretization of the har- 
monic oscillator equation can also be used to in the integration of some partial 
differential equations by Fourier transformation [6]. 

The central topic of our paper is another fruitful direction of using exact 
integrators, namely the so called locally exact discretizations [HI E] ■ 

2.2 Local exactness 

Motivated by the results of [HI IS] we propose the following definition. 

Definition 2.3. A numerical scheme Xn+i = $(a:„,/i„) for an autonomous 
equation X = F{x) is locally exact atx if its linearization around x is identical 
with the exact discretization of the differential equation linearized around x. 

The simplest choice is x = xq, where V'{xq) = (small oscillations around 
the stable equihbrium). In this case 6n does not depend on n. The resulting 
scheme, known as MOD-GR (compare [IHl [H]), was first presented in [S]. 
In [9] we considered the case x = Xn (GR-LEX) and its symmetric (time- 
reversible) modification x = |(x„ + Xn+i) (GR-SLEX). Note that in both 
cases X changes at every step. The scheme MOD-GR is locally exact at a 
stable equilibrium only, GR-LEX is locally exact at Xn (for any n), and, 
finally, GR-SLEX is locally exact at |(x„ + Xn+i) (for any n). In the case of 
implicit numerical schemes the function $(a;„,/i„) is, of course, an implicit 
function. 

Definition 2.4. A numerical scheme Xn+i = ^{Xn,hn) for an autonomous 
equation x = F{x) is locally exact if there exists a sequence Xn such that 
Xn ~ Xn = 0{hn) and the scheme is locally exact at Xn (for any n). 

Therefore GR-LEX and GR-SLEX are locally exact. Similar concept 
("linearization-preserving" schemes) has been recently formulated in p3] 
(Definition 2.1). An integrator is said to be linearization-preserving if it 
is linearization preserving at all fixed points. This definition is weaker than 
our Definition 12.41 All locally exact schemes are linearization-preserving but, 
in general, linearization-preserving schemes do not have to be locally exact 
in our sense. For instance, the scheme MOD-GR (see [8l[T0]) is linearization- 
preserving (provided that V{x) has only one stable equilibrium) but is not 
locally exact. The problem of finding the best sequence x„ for a given nu- 
merical scheme seems to be interesting but has not been considered yet. 



2.3 Exact discretization of linearized equations 

As an immediate consequence of Definition 12.41 we have that the exact dis- 
cretization of the hnearization of a given nonhnear system is locally exact. 
We confine ourselves to autonomous systems of the form 

X = F{x) (2.5) 

where x{t) G W^. If x is near a fixed x, then (12. 5p can be approximated by 

e = F\x)^ + F{x) (2.6) 

where F' is the Frechet derivative (Jacobi matrix) of F and 

^ = x-x. (2.7) 

The exact discretization of the approximated equation (12.61) is given by 

Ui = e'"^'^'^C^ + (e'^"^'^^) - /) {F'ix))-' Fix) , (2.8) 

provided that det F'{x) ^ 0. This assumption, however, is not necessary 
because function 

ipiix) := ^—^ (and cp^iO) := 1) (2.9) 

X 

is analytic also at x = 0. 

Proposition 2.5. The exact discretization of the linearization of x = F{x) 
around x is given by: 

Xn+i -x = e'^"^'^^) [xn -x) + h^ip,{h^F' (x)) F{x) . (2.10) 

The scheme Ii2.10\) is locally exact. 

If we put X = Xn into (I2.1QP (thus changing x at each step), then ^^ = 
and we obtain: 

x„+i =Xn+ (e'^"^'(-") - /) {F\xn)r'F{Xn) (2.11) 

This method is well known as the exponential difference equation (see 
formula (3.4)) or, more recently, as the exponential Euler method 



3 Locally exact modifications of popular one- 
step numerical schemes 

In order to illustrate the concept of locally exact modifications, we apply this 
procedure to the following class of numerical methods: 

Vn+l-Vn = h^iVn^yn+l) , (3-1) 

where \E'(y,t/) = F{y) and y = F{y) is equation to be solved. This class 
contains, among others, the following numerical schemes: 

• exphcit Euler scheme (EEU): ^(|/n)l/n+i) = FiVn) , 

• imphcit Euler scheme (lEU): ^(|/„,|/„+i) = F(y„+i) , 

• imphcit midpoint rule (IMP): ^(t/„,t/„+i) = F( ^"+^"+^ ) , 

• trapezoidal rule (TR): *(2/„,t/„+i) = i (F(y„) + F(y„+i)) . 

A natural non-standard modification of (13 .ip can be obtained by replacing 
/i by a matrix 5. Thus we consider the following familiy of modified numerical 
schemes: 

l/„+i-!/„ = '5(2/„,/in)^(2/„,l/„+i) , (3.2) 

where y„ is described in Definition 12.41 and S{yn, hn) is a matrix- valued func- 
tion. We assume consistency conditions: 

lini^i^^=/, M/(y,t/) = F(t/), ^,{y.,y)+^^{y,y) = F\y), (3.3) 

where J is m x m identity matrix, \E'i, \E'2 are partial Frechet derivative of ^ 
with respect to the first and second vector variable, respectively (thus \l'i, "^2 
are matrices). We also denote 

^ = ^iyn,yn) , ^1 = ^l(t/n, t/n) , ^2 = ^2(j/n,t/n) . (3.4) 

Theorem 3.1. Numerical scheme li3.S\) . where yn is any sequence described 
in Definition 2.4 and "^ satisfies ( 15*. 3\} . is locally exact for 



S{yn, hn) = hnipiihuF'iyn)) (/ + kn^ 2MhnF' (yn))) ^ . (3.5) 



Proof: The exact discretization of the Unearization of equation y = F[y) is given 
by (|2.8p . The hnearization of the scheme (j3.2p (at y„) reads 

Vn+l -Vn= 5{^iVn + ^2Vn+l + ^) , (3.6) 



where I'n = Vn —Vn, S = S{yn, hn) and we use ()3.4p . Identifying ()3.6p with ()2.8p 
and assuming invertibihty of / — 5^2 we get a system of two equations: 

(3.7) 
{I-6^2r^S^ = hnMhnF')F , 

where F = F{yn) and F' = F'{yn)- Equations (|3.7p imply 

8 i^i + ^ae'^"^') = e^-^' - I , 
^ ^ (3.8) 

(J ('I' + /l„*2¥'l(/ini^')^) = Kipi{hnF')F . 

Using (fOj) and ([Q]) (i.e., ^ = F, ^^ + ^2 = F'), and replacing e''"'^' by 
/ + hnF'ipi{hnF'), we rewrite (|3.8p as follows 

S{I + hn^2MhnF')) F' = K^i{KF')F' , 

(3.9) 
5{I + hn^2Vl{hnF')) F = K^i{KF')F . 

One can easily see that 5 given by p.Sp satisfies simultaneously both equations 
(j3.9p (actually (j3.5p is also necessary provided that F' is invertible). □ 



By straightforward calculation we can express matrix 5 in the following 
(equivalent) form: 

h F' / 1 - - h F'\^^ 

<J(^n,/in) = Manhc^(/ + -/i„(^2-^i)tanhc^j , (3.10) 

where tanhc(2;) = z'^ tanh(2;) is analytic at z = 0. For small /i„ we have 

S{Vn, K) = hj + \hl{^^ - ^2) + 0{hl) . 

Therefore 5 exists for sufficiently small hn- In other words, numerical scheme 
(13.21) admits a locally exact modification for sufficiently small /i„. 



Locally exact explicit Euler scheme 

Specializing "^{ymVn+i) = ^iVn) we obtain a class of locally exact modifica- 
tions of the explicit Euler scheme 

a:„+i = x„ + (e'^"^'^^) - l){F'{x))-'F{x^) . (3.11) 

In this case the choice x = Xn seems to be most natural. 

Proposition 3.2. Locally exact modification of the explicit Euler scheme 
(EEU-LEX) IS given by 

a:„+i = rr„ + [e^-^'^--) - l) {F'{x^))-'F{x^) . (3.12) 

Note that (KT2^ coincides with (KTl} ). 

Locally exact implicit Euler schemes 

For ^(2/„,2/„+i) = F{yn+i) we get 

Xn+i =Xn+{l- e-'^"^'(^)) {F\x))-^F{Xn+i) , (3.13) 

or, in an equivalent way, 

a;„+i =Xn + hn ipi{-hnF'{x)) F{Xn+l) ■ (3.14) 

In this case it is not clear what is the most natural identification. We can 
choose either x = Xn, or x = Xn+i- 

Proposition 3.3. Locally exact modification of the implicit Euler scheme is 
given either by 

Xn+i =Xn+{l- e-'^"^'(^")) {F'iXn)r'F{Xn+i) , (3.15) 

which will be called IE U- LEX, or by 

Xn+i = x„ + (/ - e-'^-^'(=^"+^))(F'(a:„+i))-iF(a:„+i) , (3.16) 

called lEU-ILEX. 

Both locally exact implicit Euler schemes are of second order. lEU-ILEX 
is a little bit more accurate, but lEU-LEX is more effective when computa- 
tional costs are taken into account (see Section [7]). 
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Locally exact implicit midpoint rule 

We take ^(2/„,|/n+i) = F( ^"+^"+^ ). Then, using f l310|) . we get 



.u.^^-i..^ulhnF'ix) 



<5„ = 2{F'{x))-Hanh ( ^^^^^ 1 • (3.17) 

Thus locally exact modification of the implicit midpoint rule is given by 

A f^ , hnF'{x)\ fXn+l+Xn\ .„ ^ „. 

x„+i - a;„ = /i„ I tanhc 1 F\ 1 (3.18) 

The most natural choice oix seems to be at the midpoint, x = ^ (x,„ + x„+i). 
The obtained scheme will be called IMP-SLEX. However, in order to diminish 
the computation cost, the choice x = Xn, i.e., IMP-LEX, can also be consid- 
ered (because then Jacobian F'{x) is evaluated outside iteration loops). 

Locally exact trapezoidal rule 

Takin g into account ^(t/„,t/„_^.J = | (F(t/„) + F(y„+i)) and p.iup, we get 
(I3.17p . as well. Thus locally exact modification of the trapezoidal rule reads 



In i tanhc^^-- j 



x„+i -Xn = K[ tanhc ^^ ] — '- ^ — - . (3.19) 



There are two natural choices of x. Either (in order to minimize the compu- 
tational costs) we can take x = Xn, obtaining TR-LEX, or (in order to get 
time-reversible scheme TR-SLEX) we can take x = ^ (a;„ +Xn+i)- 

4 Locally exact discrete gradient schemes for 
multidimensional Hamiltonian systems 

In this section we construct energy-preserving locally exact discrete gradient 
schemes for arbitrary multidimensional Hamiltonian systems in canonical 
coordinates: 

.;. _ dH .^ _ dH 

where k = 1, . . . ,m. We obtain two different numerical schemes using either 
symmetric discrete gradient or coordinate increment discrete gradient. 



4.1 Linearization of Hamiltonian systems 

We denote x = {x^, . . . ,x"^)'^, p = {p^, ■ ■ ■ ,p"^), etc. The linearization of 
(14. ip around x,p is given by: 

^ = Hp + Hpx^ + HppTi , 

(4.2) 
V = —Hx — Hxxi — HxpT] , 

where ^ = x — x, v) = p — p, Hp is a vector with components |^, Hpx is 

f:j2 Tj 

a matrix with coefficients g jQ^k , etc., and derivatives of H are evaluated at 
x,p. Note that m x m matrices H^x, Hxp, Hpx, Hpp satisfy 



PP ~ PP ' XX ~ ^XX ) ^xp ~ ^PX l4.dj 



Equations ( 14. 2 p can be rewritten also as 

where 

F = ( ^^ ] , F' = ( ^^^ ^^ ] . (4.5) 

\^ -^x J \ ^xx -^xp J 

Corollary 4.1. The exact discretization of linearized Hamiltonian equations 
lli4-4^ is given by 

Proof: The exact discretization of (j4.4p is given by ()4.6p which follows immediately 
from Proposition 12.51 □ 



4.2 Discrete gradients in the multidimensional case 

Considering multidimensional Hamiltonian systems we will denote 
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where x E M"^, p E M"*, y E W"^, etc. In other words, 

A discrete gradient, denoted by Vi/(t/„,t/„+i), 

- ( AH AH AH\_(AH AH\ 

\A^'Af''"'A^)^{'A^'^) ^ ' 

is defined as an R^™- valued function of t/„,y„+i such that [T^ [22]: 

k=i ^ (4.9) 

VH{y,y) = {H,,Hp) , 

where ifi, -f^p are evaluated at y = (a;,p) and we define 

VH{y,y) = \imVH{y,y) (4.10) 

when necessary. It is well known (see, e.g., [121 1131 1201 [22]) that any discrete 
gradient defines an energy-preserving numerical scheme 

yn+i = yn + hVH . (4.11) 

Discrete gradients are non-unique, compare [T21 [T^l [22]. Here we con- 
fine ourselves to the simplest discrete gradients. We consider coordinate 
increment discrete gradient, first proposed by Itoh and Abe [IS], and its 
symmetrization (see below). The coordinate increment discrete gradient is 
defined by: 



AH 


H{yl_^^,ylyl... 


,y?r)-H{ylylyl--- 


, yin 


Ay' 
AH 


H{yl^^,yl^^,yl. 


yl+i - yl 
..,y^)-i/(yi+i,2/2,.. 


7 


Ay^ 




yl+i - yl 


•) 


AH 


_ // (2/^+1,2/^+1,.., 


.,y^Ti) -^fci, 2/^+1, • 


•-,2/^) 



(4.12) 



At, 2m ,,2m ,,2m 

^2/ 2/n+l 2/n 
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where, to fix our attention, y is defined by f l4.7p . The corresponding numerical 
scheme fl4.1ip will be called GR-IA. In fact, we may identify with y any 
permutation of 2m components x^ ^jP . Thus we have (2m)! discrete gradients 
of this type (some of them may happen to be identical). 

Having any discrete gradient we can easily obtain the related symmetric 
discrete gradient 



VsH{yn,yn+i) = - (V//(y„,|/„+i) + Vi/(y„+i,y„)) . 



(4.13) 



One can easily verify that VgH satisfies conditions (14. 9 p provided that they 
are satisfied by VH. Discrete gradient scheme obtained by symmetrization 
dUSD from GR-IA will be called GR-SYM. 



4.3 Linearization of discrete gradients 

In order to construct locally exact modifications we need to linearize discrete 
gradients defined by (I4.12p and (I4.13p . 



Lemma 4.2. Linearization of the coordinate increment discrete gradient 
Ijl4-i2 ) around y yields 



1 



Hy + - {AUn+l + BUn) 



VH{yn,yn+i] 
where we denoted i/„ = Vn — V, (md 



(4.14) 



A 



Hy2yl 





-H 2 2 

2 y y 








liy2m—lyl 



\Hy 



T2y2m — ly2 
Ily2my2 








/ ^Hylyl Hyly2 



B 







-H 2 2 

2 y y 












~ liy2m—ly2m — l 
11 y2m y2m—l 



\ 



T: riy2iny2m 



llyly2m—l 
rly2y2m-l 



llyly2m 
Ily2y2m 



7^ ^n. y2m — ly2m — l 



lly2m — l y2m 

"^ Ily2my2m 



(4.15) 



where H, 



y]yk 



dyidy^ 
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Proof: We denote yi = (y^+i, . . . , y^+i, yn''"\ • • • , V^TV ■ In particular, y° = y„ 
and y^™ = 2/„,+i. Expanding H{yn) around yi-i~ , we get 

Hiiii) = F(yr')+^,.(yr')fci-y^)+^i:^,.,.(yr')(y^+i-y'J'+--- (4.i6) 

Hence 

Then, from the definition of i^n, we have 



J j J J 

2/n+l ~ Vn — ^n+1 ~ ^n ) 



,V-1 - ,-7 _L h.l . ,.^-1 ,.i z,2m 



^n — 2/ + ( '^n+1) • • • ) ^n+1' '^n' • • • ' "n / ' 



(4.18) 



and, taking it into account, we rewrite (j4.17p as 

Ai? ^^^ ^"^ 1 • • 

^"7 = Hyj (y) + Y, Hyjyk (y) i^n+i + Yl ^y'y" (^) ""n + ^Hy.yj (y) {u^^^ -i^i) + ..- 

^ k=l k=j 

(4.19) 
which is equivalent to (I4.14P , (|4.15p . □ 

Lemma 4.3. Linearization of the symmetric discrete gradient yields 

VsH{yn, ?/„+i) ^ Hy + -Hyy (z/„ + i/„+i) , (4.20) 

where Hyy is the Hessian matrix of H , evaluated at y. 

Proof: We observe that A + B = Hyy and then we use (|4.13p and (|4.14p . □ 

4.4 Conservative properties of modified discrete gra- 
dients 

In the one-dimensional case locally exact modifications are clearly energy- 
preserving, compare [9], [10]. In the general case, conservative properties are 
less obvious. In this section we present several useful results. 
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Lemma 4.4. We assume that a 2m x 2m matrix A (depending on h and, 
possibly, on other variables) is skew- symmetric (i.e., A'^ = —A) and 

Then, the numerical scheme 

Vn+l -Vn = AVif , (4.22) 

where yn is defined by ( [^.Tp and VH satisfies ^.9^ , is a consistent integrator 
for ([^. j|j preserving the energy integral up to round-off' error. 



Proof: The consistency follows immediately form (|4.2ip . The energy preservation 
can be shown in the standard way. Using the standard scalar product in R^™, we 
multiply both sides of (lO^ll by VH 

{VH I y„+i - y„) = {VH \ AVH) . (4.23) 

By virtue of (14. 9 p the left-hand side equals H{yn+i) — H{yn). The right hand side 
vanishes due to the skew symmetry of A. Hence -ff(y„+i) = H{yn). □ 

Remark 4.5. Scheme i4.22^ for A = hS becomes a standard discrete gradi- 
ent scheme. 

Lemma 4.6. We assume that 2m x 2m matrix 6 is of the following form 

5 —O \ T T ■• ^ 



(where 5, o , p are m x m matrices) and VH is (any) discrete gradient. 
Then, the numerical scheme given by 

|/„+i -yn = eSVH (4.25) 

preserves the energy integral exactly, i.e., H{Xn+i,yn+i) = H{Xn,yn)- 

Proof: We observe that dS is skew-symmetric, and then we use Lemma |4.4[ □ 

We easily see that f l4.25p reduces to the standard discrete gradient scheme 
when we take 6 = hnl (i.e., 6 is proportional to the unit matrix). 
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Lemma 4.7. // 6 is of the form ( [^ . 24\ ) and z ^^ f{z) is any analytic 
function, then f{9) is also of the form ( [^J 



Proof: First, we will show that the conditions (j4.24p are equivalent to 

d'^ = S-^dS . (4.26) 

Indeed, assuming a general form oi 6, e.g., 6 = i | we see that the con- 

straint ()4.26p is equivalent to f = 6'^ , p^ = —p and <t^ = —a. Then the proof is 

oo 

straightforward. Assuming f{z) = y^OnZ^, we obtain 

k=l 

/ oo \ 



^ a^e'' =Y,ak {S-'dS) ' = S-UY, ake" S , (4.27) 

\A:=1 / fc=l \fc=l / 

i.e., fief = s-^f{e)s. n 



Corollary 4.8. // 6'^ = S ^6S and f is an analytic function, then the 
scheme t/n+i — Vn = f{G)SVH exactly preserves the energy integral H . 

4.5 Locally exact symmetric discrete gradient scheme 

We begin with the symmetric case because the coordinate increment dis- 
crete gradient case is more difficult. In the symmetric case a locally exact 
modification is derived similarly as in the one-dimensional case. 

Proposition 4.9. The following modification of a symmetric discrete gradi- 
ent scheme is locally exact at y: 

Vn+l - Vn = OnSW sH , (4.28) 

where 

9n = 2{F')-Hanh^ , (4.29) 



an 



d F' , given by ( [^.5] j, is evaluated at y = y. 
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Proof: We are going to derive (j4.29p . assuming that 9n depends on y and h. By 
virtue of Lemma 14.31 the Hnearization of (j4.28p is given by 

l^n+l -l'n= dnS i Hy + -Hyy{Vn+l + Vn)\ • (4.30) 

Taking into account (j4.5p we transform (|4.30p into 

(l - \6nF'\ Vn+1 = f / + \6nF'\ V^ + 0nF . (4.31) 

The scheme (|4.28p is locahy exact if and only if (j4.3ip coincides with (j4.6p . There- 
fore, we require that 



(/ - ie„F') e'"^' = I + \dnF' 
I - ^^dnF^ (e^-^' - l) {F'r'F = dnF . 



(4.32) 



(4.33) 



From equation ()4.32p we can compute On which yields (j4.29p . Equation (j4.33p is 
automatically satisfied provided that (14.32P holds. □ 

Taking the syinmetrization of the Itoh-Abe discrete gradient (compare 
(I4.13P ). modifying it according to (14.281) and choosing y = ?/„ we get a scheme 
called GR-SYM-LEX, while for t/ = i (y„ + y^+i) we obtain GR-SYM-SLEX. 



Proposition 4.10. The numerical scheme lli4-28[ ) with On given hy ^4-29^ is 
energy-preserving. 

Proof: We have F' = SHyy. Therefore, (F')^ = —HyyS = —S^^F'S, and, as a 
consequence 

{{F'ff = S-\F'fS . (4.34) 

It means that {F')'^ has the form (j4.24p . The formula (14.291) expresses dn as an 
analytic function of {F')"^. Finally, we use Lemma |4.7[ □ 

In the one-dimensional case {F'Y' is proportional to the unit matrix which 
essentially simphfies arguments presented in this section, see [91 ITU]. 
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4.6 Locally exact coordinate increment discrete gradi- 
ent scheme 

The symmetric form of the discrete gradient leads to a simple form of the 
locally exact modification. It turns out, however, that starting from the 
simplest form of the discrete gradient, namely coordinate increment discrete 
gradient GR-IA, we also succeed in deriving corresponding locally exact mod- 
ifications. 

Proposition 4.11. The following modification of the coordinate increment 
discrete gradient scheme is locally exact at y: 



l/n+l -Vn = OnSVH , 

where VH is given by ^.IS^ ), 

2(sR + F'coth^ 



( 



R 



-H 1 2 



H 2 1 



— Ily2y2m-1 



~ 12 yl y2m 
~ lly2y2Tn 



\ 



Hy2m -lyl Hy2m -\y1 

\ IIy2myl IIy2my2 



' £1 y2m — l y2m 



Ily2my2m — 1 



F' and R are evaluated at y = y. 



(4.35) 



(4.36) 



F' is given by ( [^.5| j (i.e., F' = SHyy), and, finally R = A — B, i.e.. 



(4.37) 



Proof: We are going to derive ()4.36p . assuming that 9n depends on y and h. By 
virtue of Lemma 14.21 the linearization of (j4.35p is given by 



l^n+l — J^n = OnS{Al>n+l + BVn) + OnSHy 

Hence, taking into account that SHy = F, 

(/ - OnSA) Vn+1 = {I + 6nSB) V^ + BnF . 



(4.38) 



(4.39) 



The scheme (j4.35p is locahy exact if and only if (j4.39p coincides with (j4.6p . There- 
fore, we require that 



{I-dnSA)e''"^ =I + dnSB , 



(4.40) 
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(I - OnSA) [e^"^' - Ij {F'r^F = OnF . (4.41) 

Inserting (I4.40p into (14.41 1) we get 

dnSiB + A){F'r^F = dnF , (4.42) 

which is identically satisfied by virtue oi A + B = Hyy = S^^F', compare (j4.15p . 
The remaining equation, (j4.40p . defines 0„: 



er 



(sAe''"^' + Sb) = e''"^' - I . (4.43) 



In order to simplify ()4.43p we introduce R = A — B {R is antisymmetric because 
B = A^^). Then, taking into account SA + SB = F', we get 

SA = -F' + -SR , SB = -F' - -SR . (4.44) 

2 2' 22 ^ ^ 

Substituting it into (j4.43p we complete the proof. □ 

Choosing y = Vn ^e obtain GR-IA-LEX, while for y = | (t/„ + Vn+i) we 
get GR-IA-SLEX. 



Proposition 4.12. The numerical scheme { 4-35 ) with On given by { 4-36 ) is 
energy-preserving, i.e., H{yn+i) = H{yn). 



Proof: We have 



T\ -1 



el = 2UsRf+^F'coth^^ j =s-'enS 



(4.45) 



because 

{SRf = {-R){-S) = S-^ {SR) S (4.46) 

and, by virtue of Lemma 14.71 

F'coth^] = S-^ (f' coth^^ S , (4.47) 

where we took into account (j4.34p . Then, we use Corollarv 14.81 □ 



Locally exact modification fl4.35p for Hamiltonian H = T{p) + V{x) with 
two degrees of freedom is determined by 



dn = /in.tanhc 



K.F' f 1 hnF' 

— — I / H — hnS Rtanhc— — 

2 V 2 2 



-1 



(4.48) 



where 



R 



( -y,i2 \ 

l/,i2 

-T,i2 

V T,i2 / 



(4.49) 



In numerical experiments we use T{;p) = |p^ and a potential depending 

only on r = |a;|. In this case T,i2 = and l^,i2 = ^^ il^yr) rr- After 
straightforward calculations we reduce (I4.35P to: 



^71+1 '^n — riiy 



Pn+1 + Pn 



Pn+1 - Pn = -hD VV - h^det D) V,i2 ( \ J 



where Q"^ = V,xx and D = tanhc^. 



(4.50) 



5 Order of considered numerical schemes 

First of all, we recall the Taylor expansion of the exact solution of equation 
X = F{x) (compare, e.g., [21 1 



x{t + h)=x + hF+ h^F'F + ^h' {F"iF, F) + {F'fF) + ... (5.1) 

where all terms on the right-hand side are evaluated at t. Note that F is a 
vector, F' is a matrix, F" is a vector-valued bilinear form, etc. Then, x{t) and 
x(t+h) will be identified with Xn and Xn+i, respectively. To simplify notation, 
in this section we replace hn by h. Computing the order of numerical schemes 
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(5.2) 



we use expansions 

F{xn+i) = F + hF'h + h^ (F'b2 + liF"ih, b,)] 

+h' (F'bs + F"ibi, 62) + liF'"ibi, 61, 61)") + . . . 

F (^^^i±^^ =F + \hF\ + h^ {^-F% + 1f"(6„ 6,)) 

+h' Qf'63 + lF"ib„ b,) + ^F'"{b,, 61, 61)) + . . . 

where bk are coefficients of the Taylor expansion of ^{Xn, h), i.e., 

Xn+i = Xn + bih + b2h^ + b^h? + . . . (5.3) 

In the case of locally exact schemes we use also the following Taylor expan- 
sions 

f,(hF) = / + \hF' + j,h\Fr- + ^h\F'f + ... 

substituting (15. 2 p if F' is evaluated at x„+i or | (x„ +Xn+i)- The final re- 
sult of this analysis (expansions (15. 3p for particular numerical schemes) is 
presented below. 

First order schemes 

• Explicit (forward) Euler scheme, EEU, 

Xn+l =Xn + hF . (5.5) 

• Implicit (backward) Euler scheme, lEU, 

x„+i =Xr, + hF + h^F'F + ... (5.6) 
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Second order schemes 

• Locally exact explicit Euler scheme {x = x„), EEU-LEX, 

x„+i =Xn + hF+ h^F'F + Ih'iF'fF + . . . (5.7) 

• Locally exact implicit Euler scheme {x = x„), lEU-LEX, 

a:„+i =Xn + hF + h'F'F + h' (^{F'fF +^F"{F, F)] + . . . (5.8) 

• Locally exact implicit Euler scheme {x = Xn+i), lEU-ILEX, 

Xn+i =Xn + hF + hi^F'F + h? {^-{F'fF + 1f"(F, F)") + . . . (5.9) 

• Implicit midpoint rule, IMP, 

Xn+i=Xn + hF + hi^F'F + h'{hF'fF+^-F'\F,F)\+... (5.10) 

• Locally exact implicit midpoint rules IMP-LEX, IMP-SLEX (/i-expan- 
sions for x = x„ and x = |(x„ + Xn+i) are identical up to the third 
order) : 

Xn+i = Xr. + hF + ^h'^F'F + h^ (liF'fF + ^F"{F, F)\ + . . . (5.11) 

• Trapezoidal rule, TR, 

Xn+i = Xn + hF + ^^F'F + h^ (liF'fF + ^F"{F, F)\ + . . . (5.12) 

• Locally exact trapezoidal rules TR-LEX, TR-SLEX (/i-expansions for 
X = ^{Xn + Xn+i) and X = Xn are identical up to the third order): 

Xn^i = Xn + hF + h'^F'F + h' (^{F'^F +^F"{F, F)\ + . . . (5.13) 
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Gradient schemes 

Discussing gradient schemes we confine ourselves to Hamiltonians of tlie form 
H = T{p) + V{x). First, we expand a discrete gradient of V with respect to 
Ax-' (jth component of a;„+i — x„, compare (14. 8p ): 

^ = V,,, +^A,kAx' + ^B,^,Ax''Ax'' + ^C,kf.uAx''Ax^Ax'' + . . . (5.14) 

where Aj^, Bj^^, and Cjk^u are x„-dependent coefficients. In particular: 

• Itoh-Abe gradient , GR-IA, 

Ajj = V,^j^j, Ajk = ij<k), Ajk = 2V,^,^k {j>k), (5.15) 

• any symmetric discrete gradient, compare (I4.13p . 

Ajk = y.xix^ , (5.16) 

• symmetrization of the Itoh-Abe gradient, one degree of freedom, 

^11 = y^xx , Bill = V,xxx , Ciiu = V,xxxx , (5-17) 

• symmetrization of the Itoh-Abe gradient (GR-SYM), two degrees of 
freedom, 



Bill — V,xxx , Bii2 — -V,xxy , Bi22 — T^^ixyy ■, 

^ ^ (5.18) 

3 3 

B222 = ^lyyv 5 B212 = -V,xyy , -B211 = —V,xxy ■ 

Having (I5.14p we substitute expansion of Ax^ with respect to h in the 
considered discrete gradient scheme. As a results we obtain Taylor series for 
Xn+i and comparing them with (15. ip we arrive at the following conclusions 
(assuming the generic case, because for some exceptional V, e.g., V linear or 
quadratic in x^ , the order can be higher). 

First order schemes 

• Discrete gradient schemes such that Aj^ 7^ Vj^jx'', in particular: the 
Itoh-Abe scheme GR-IA. 
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Second order schemes 

• Discrete gradient schemes such that Ajk = V,x3x''i in particular: sym- 
metric gradient schemes (including GR-SYM) and one-dimensional dis- 
crete gradient method. 

• Locally exact {x = x„, i.e., LEX) modifications of sjTiimetric discrete 
gradient schemes such that Bj^^, ^ V^^jxt^x"- In particular, locally ex- 
act modification of symmetrization of the Itoh-Abe discrete gradient 
scheme (GR-SYM-LEX), multidimensional case (i.e., two degrees of 
freedom, at least). 

• Locally exact time reversible {x = ^ {Xn +Xn+i), i.e., SLEX) modifica- 
tions of symmetric discrete gradient schemes such that -Bj^^ 7^ V^x^xi^x''- 
In particular, locally exact time reversible modification of symmetriza- 
tion of the Itoh-Abe discrete gradient scheme (GR-SYM-SLEX), mul- 
tidimensional case. 

• Locally exact {x = aj„) modification of the Itoh-Abe discrete gradient 
scheme (GR-IA-LEX). 

• Locally exact time reversible (S = | (x„ +a;„+i)) modification of the 
Itoh-Abe discrete gradient scheme (GR-IA-SLEX). 

Third order schemes 

• Locally exact (LEX: x = x„) modifications of symmetric discrete gra- 
dient schemes such that Bj^^, = V^xixt^x"- In particular, GR-LEX in 
one- dimensional case, see [TU] . 

• Locally exact time reversible (^ = | (a;„ +Xn+i), i.e., SLEX) modifica- 
tions of symmetric discrete gradient schemes such that Bj^^, = V^^j^i^xu 

and iyjknu 7^ ^ ixix'^xt^x^ • 

Fourth order schemes 

• Locally exact time reversible (SLEX: x = \ (x„ +a;„+i)) modifications 
of symmetric discrete gradient schemes such that Bj^^ = V^xixt^x" and 
CjkiJLv = ^jxJx'^xfa"- In particular, GR-SLEX in one-dimensional case, 
see [in]. 
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6 Linear stability 

Locally exact numerical schemes have excellent qualitative behaviour around 
all fixed points of the considered system. 

Proposition 6.1. // the equation x = F{x) has a fixed point at x = x, then 
all its locally exact discretizations have a fixed point at Xn = x, as well. 

Proof: If F{x) = 0, then equation (|2.10p becomes 

Xn+l-X = e^-^'^-^Hxn-x) . (6.1) 

We require that the scheme Xn+i = ^{Xn,h) is a locaUy exact discretization of 
F{x) = 0, i.e., its hnearization, given by 

Xn+i = $(x, h) + ^'{x, h){Xn - x) , (6.2) 

coincides with (16.11). Hence 



e 



hr^F'(x) 



<^'{x,h) , $(x,/i) =x , (6.3) 



which means that a; is a fixed point of the system Xn+i = ^{xn, h). D 

Stability of numerical integrators can be roughly defined as follows: "the 
numerical solution provided by a stable numerical integrator does not tend to 
infinity when the exact solution is bounded" (see [3], p. 358). The integrator 
which is stable when applied to linear equations is said to be linearly stable. 
We may use the notion of A-stability (see, e.g., [H]): an integrator is said to 
be A-stable, if discretizations of all stable linear equations are stable as well. 

Making one more assumption (quite natural, in fact) that the discretiza- 
tion of a linear system is linear, we see that locally exact integrators are 
linearly stable. Indeed, solutions of any locally exact discretization have the 
same trajectories as corresponding exact solutions. 

Corollary 6.2. Any locally exact numerical scheme is linearly stable and, in 
particular, A-stable. 

What is more, a locally exact discretization yields the best (exact) simula- 
tion of a linear equation in the neighbourhood of a fixed point. In particular, 
locally exact discretizations preserve any qualitative features of trajectories 
of linear equations (up to round-off errors, of course). 
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Numerical experiments show that locally exact schemes are exceptionally 
stable also for some simple nonlinear systems [51 13 dUI, but we have no 
theoretical results concerning the stability (e.g., algebraic stability [H]) in 
the nonlinear case. 

In order to illustrate these general results we will apply four schemes 
presented above to a linear equation x = Ax. Locally exact explicit Euler 
scheme and locally exact implicit Euler scheme yield, respectively, 

Xn+i -Xn= (e^"^ - /) x„ , (6.4) 



Xn+l -Xn= {l-e ^"^ - /) Xn+1 ■ (6.5) 

Implicit midpoint and trapezoidal rules yield an identical equation, namely 

Xn+i -Xn= i tanh ^- J {x^+i + Xn) ■ (6.6) 

Simple calculations show that all resulting equations reduce to the exact 
discretization: Xn+i = e'^"^Xn- In particular, if real part of any eigenvalue of 
A is negative, then x„ — )■ for n — )• oo. A-stability is evident. 



7 Numerical experiments 

In order to illustrate performance of locally exact numerical schemes we con- 
sider circular orbits for the Hamiltonian 

H(x,p) = -\p\^ + -\x\^-—\xf, (7.1) 

where the exact solution can be easily found. Namely, initial conditions 
Xq = (2;o;0), po = (0,Xoa/1 — O.lxo) yield a circular orbit of radius R = Xq. 
Circular orbits exist for i? < 10. We made numerical tests for R= 0.2, R = 1 
and R = 5 (coresponding periods of exact solutions are given by: 6.347, 6.623 
and 8.886, respectively). While solving implicit algebraic equation for Xn+i 
we used the fixed point method. Iterations were made until the accuracy 
3 ■ 10~^® was reached. The maximal number of iterations was limited to 20. 
All figures present global error at t = 12.5. 

Locally exact modifications improve the accuracy but are quite expensive. 
In our numerical experiments we use different time steps for different schemes 
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to assure the same computational cost (i.e., at all figures the computational 
cost in every column is the same). The numerical cost is estimated as a 
number of function evaluations. In other words, we use the scaled time step 
h such that h = Xh, where A depend on numerical scheme and (to some 
extent) also on h. In order to evaluate the parameter A we start from A = 1, 
estimate the computational costs of considered schemes and multiply h by the 
relative computational cost. The procedure is repeated until computational 
costs become equal with accuracy about 1%. Approximate values of the 
parameter A are given in figure captions. 

In the case of Euler schemes (Figs. [1] and [2]) the advantage of locally exact 
modifications is obvious. Although the cost of locally exact modifications is 
considerably higher, they perform much better than both standard Euler 
schemes. Actually the accuracy of all 3 modifications (taken with the same 
time step) is similar. However, taking into account the computational cost, 
we see that the scheme EEU-LEX is the best modification. 

Locally exact modifications of implicit midpoint and trapezoidal rules 
increase have higher accuracy only for orbits with a relatively small radius 
(Figs. OandH]). The best results are produced by scheme IMP-LEX. For 
larger orbits the unmodified implicit midpoint rule is most accurate. Similar 
situation has place for gradient schemes. Locally exact modifications give a 
considerable improvement, by about 1-2 orders of magnitude, for trajectories 
in a quite large neighborhood of the stable equilibrium, see Figs. E] and Ei 
For large R (e.g., R = 5) locally exact modifications improve considerably 
only GR-IA. Their influence on the GR-SYM is neutral or even negative. In 
fact GR-SYM has similar accuracy for all considered orbits while its locally 
exact modifications are very accurate only for smaller values of R. 

8 Concluding remarks 

We presented a new class of numerical schemes characterized by the so called 
local exactness. This notion is known (under different names) since almost 
fifty years. The original application, see [27] , has been confined to the exact 
discretization of linearized equations, compare Section [231 We obtain in this 
way a particular locally exact integrator which has some advantages (e.g., it 
can serve as a good predictor, see [9], Section V.A). 

Our approach has two new features. First, we modify known numerical 
schemes in a locally exact way. Any numerical scheme admits at least one 
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(usually more) natural locally exact modification, see Section [31 Second, we 
try to preserve geometric properties of the original numerical scheme. This 
task is not trivial. In this paper we present one successful application: lo- 
cally exact modifications of discrete gradient methods for canonical Hamilton 
equations. We are able to exactly preserve the energy integral and, in the 
same time, increase the accuracy by many orders of magnitude. Another 
advantage is a variable time step. Unlike symplectic methods (which work 
mostly for the constant time step) discrete gradient methods admit con- 
servative modifications with variable time step. Therefore, one may easily 
implement any variable step method in order to obtain further improvement. 

In the one-dimensional case proposed modifications, although a little bit 
more expensive, turns out to be more accurate even by 8 orders of magni- 
tude in comparison to the standard discrete gradient scheme, see [9l |T0]. In 
multidimensional cases the relative cost of our algorithm is higher, but still 
our method is of great advantage for orbits in the neighborhood of the stable 
equilibrium. The accuracy is increased by 1-2 orders of magnitude. We point 
out that in most cases locally exact modifications do not change orders of 
modified schemes. 

We presented locally exact modifications from a unified theoretical per- 
spective. There are many possible further developments. First of all, we plan 
to apply our approach to chosen multidimensional problems, testing the accu- 
racy of locally exact schemes by numerical experiments. Then, we would like 
to extend the range of applications. It would be natural to consider locally 
exact modifications of Runge-Kutta schemes and G-symplectic integrators 
[2]. However, our first attempts seem to suggest that this is a challanging 
problem. Modifications proposed in this paper contain exponentials of vari- 
able matrices. Similar time-consuming evaluations are characteristic for all 
exponential integrators and in this context effective methods of computing 
matrix exponentials have been recently developed [151 126] . 

Throughout this paper we assumed the autonomous case, x = F{x). 
The extension on the non-autonomous case can be done along lines indi- 
cated already in Pope's paper [27] . A separate problem is to obtain in this 
case any locally exact modification with geometric properties. Another open 
problem is the construction of locally exact (or, at least, linearization pre- 
serving) defomations of generalized discrete gradient algoritms preserving all 
first integrals (see [22]). Linearization-preserving integrators form an im- 
portant subclass of locally exact schemes. It would be worthwhile to study 
linearization-preserving modifications of geometric numerical integrators, es- 
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pecially in those cases when locally exact are difficult or impossible to con- 
struct. 
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References 

[1] R.P.Agarwal: Difference equations and inequalities (Chapter 3), Marcel Dekker, New 

York 2000. 
[2] J. C. Butcher: "Numerical methods for ordinary differential equations", second edi- 
tion, Wiley & Sons, Chichester 2008. 
[3] E.Celledoni, D.Cohen, B.Owren: "Symmetric exponential integrators with an appli- 
cation to the cubic Schrodinger equation", Found. Comput. Math. 8 (2008) 303-317. 
[4] J.L.Cieslihski: "An orbit-preserving discretization of the classical Kepler problem", 

Phys. Lett. A 370 (2007) 8-12. 
[5] J.L.Cieslihski: "Comment on 'Conservative discretizations of the Kepler motion' ", 

J. Phys. A: Math. Theor. 43 (2010) 228001 (4pp). 
[6] J.L.Cieslihski: "On the exact discretization of the classical harmonic oscillator equa- 
tion", J. Difference Equ. Appl. 17 (2011) 1673-1694. 
[7] J.L.Cieslihski, B.Ratkiewicz: "On simulations of the classical harmonic oscillator 
equation by difference equations". Adv. Difference Eqs. 2006 (2006) 40171 (17pp). 
[8] J.L.Cieslihski, B.Ratkiewicz: "Long-time behaviour of discretizations of the simple 

pendulum equation", J. Phys. A: Math. Theor. 42 (2009) 105204 (29pp). 
[9] J.L.Cieslihski, B.Ratkiewicz: "Improving the accuracy of the discrete gradient 
method in the one-dimensional case", Phys. Rev. E 81 (2010) 016704 (6pp). 
[10] J.L.Cieslihski, B.Ratkiewicz: "Energy-preserving numerical schemes of high accu- 
racy for one-dimensional Hamiltonian systems", J. Phys. A: Math. Theor. 44 (2011) 
155206 (14pp). 
[11] J.L.Cieslihski, B.Ratkiewicz: "Discrete gradient algorithms of high-order for one- 
dimensional systems" , Comp. Phys. Comm. 183 (2012) 617-627. 
[12] O.Gonzales: "Time integration and discrete Hamiltonian systems", J. Nonl. Sci. 6 

(1996) 449-467. 
[13] D.Greenspan: "An algebraic, energy conserving formulation of classical molecular 

and Newtonian n-body interaction". Bull. Amer. Math. Soc. 79 (1973) 432-427. 
[14] E.Hairer, C.Lubich, G. Wanner: Geometric numerical integration: structure-preser- 
ving algorithms for ordinary differential equations, second edition. Springer, Berlin 
2006. 
[15] M.Hochbruck, Ch.Lubich: "On Krylov subspace approximations to the matrix ex- 
ponential operator", SIAM J. Numer. Anal. 34 (1997) 1911-1925. 
[16] M.Hochbruck, Ch.Lubich: "A Gautschi-type method for oscillatory second-order 
differential equations", Numer. Math. 83 (1999) 403-426. 



28 



[17] A.Iserles: "Insight, not just numbers", Proceedings of the 15th IMACS World 
Congress^ vol. II, ed. by A.Sydow, pp. 589-594; Wissenschaft & Technik Verlag, 
Berlin 1997. 

[18] A.Iserles: A first course in the numerical analysis of differential equations, second 
edition, Cambridge Univ. Press 2009. 

[19] T.Itoh, K.Abe: "Hamiltonian conserving discrete canonical equations based on vari- 
ational difference quotients", J. Comput. Phys. 77 (1988) 85-102. 

[20] R.A.LaBudde, D.Greenspan: "Discrete mechanics - a general treatment", J. Com- 
put. Phys. 15 (1974) 134-167. 

[21] D.J.Lawson: "Generalized Ruge-Kutta processes for stable systems with large Lips- 
chitz constants", SIAM J. Numer. Anal. 4 (1967) 372-380. 

[22] R.I.McLachlan, G.R.W.Quispel, N.Robidoux: "Geometric integration using discrete 
gradients", Phil. Trans. R. Soc. London A 357 (1999) 1021-1045. 

[23] R.I.McLachlan, G.R.W.Quispel, P.S.P.Tse: "Linearization-preserving self-adjoint 
and symplectic integrators", BIT Numer. Math. 49 (2009) 177-197. 

[24] R.E.Mickens: Nonstandard finite difference models of differential equations, World 
Scientific, Singapore 1994. 

[25] B.V.Minchev, W.M.Wright: "A review of exponetial integrators for first order semi- 
linear problems", preprint A^T7VC//Numerics/N2/2005, Trondheim 2005. 

[26] J.Niesen, W.M.Wright: "Algorithm 919: a Krylov subspacc algorithm for evaluating 
the (^-functions appearing in exponential integrators", ACM Trans. Math. Software 
38 (3) (2012), article 22. 

[27] D.A.Pope: "An exponential method of numerical integration of ordinary differential 
equations", Commun. ACM & (8) (1963) 491-493. 

[28] R.B. Potts: "Differential and difference equations". Am. Math. Monthly 89 (1982) 
402-407. 

[29] G.R.W.Quispel, G.S. Turner: "Discrete gradient methods for solving ODE's numeri- 
cally while preserving a first integral", J. Phys. A: Math. Gen. 29 (1996) L341-L349. 



29 



Figure 1: Error of numerical solutions at f = 12.5 for a circular orbit (R = 
0.2) as a function of scaled time step h (i.e., h = Xh). EEU: white discs 
(A = 1), lEU: white squares (A = 24,30,42,60, respectively), EEU-LEX: 
gray discs (A = 8), lEU-LEX: gray squares (A = 34,40,53,67), lEU-lLEX: 
black squares (A = 130,200,200,200). 
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Figure 2: Error of numerical solutions at t = 12.5 for a circular orbit {R = 5) 
as a function of scaled time step h (i.e., h = \h). EEU: white discs (A = 1), 
lEU: white squares (A = 25,31,42,60), EEU-LEX: gray discs (A = 8), lEU- 
LEX: gray squares (A = 34,40,52,67), lEU-ILEX: black squares (A = 200). 











n 








n 


■ 


1. 


n 


n 


■ 


■ 
o 


0.1 


s 


■ 

o 


o 

■ 




0.01 




■ 




• 


0.001 


■ 




• 




0.0001 


• 


• 







0.0005 



0.001 



0.002 



0.004 



31 



Figure 3: Error of numerical solutions at t = 12.5 for a circular orbit {R = 
0.2) as a function of scaled time step h (i.e., h = Xh). IMP: white discs 
(A = 1), IMP-LEX: gray discs (A = 1.3, 1.3, 1.3, 1.2), IMP-SLEX: black discs 
(A = 3.4,3.4,3.9,3.9), TR: white squares (A = 1.06), TR-LEX: gray squares 
(A = 1.9, 1.6, 1.8, 1.8), TR-SLEX: black squares (A = 3.5, 3.4, 4.0, 3.9). 
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Figure 4: Error of numerical solutions at t = 12.5 for a circular orbit {R = 
1) as a function of scaled time step h (i.e., h = Xh). IMP: white discs 
(A = 1), IMP-LEX: gray discs (A = 1.5, 1.3, 1.4, 1.2), IMP-SLEX: black discs 
(A = 3.4,3.4,3.9,3.9), TR: white squares (A = 1.06), TR-LEX: gray squares 
(A = 1.9, 1.6, 1.8, 1.8), TR-SLEX: black squares (A = 3.5, 3.4, 4.0, 3.9). 
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Figure 5: Error of numerical solutions at f = 12.5 for a circular orbit [R = 
0.2) as a function of scaled time step h (i.e., h = Xh). GR-IA: white discs, 
GR-IA-LEX: gray discs, GR-IA-SLEX: black discs, GR-SYM: white squares, 
GR-SYM-LEX: gray squares, GR-SYM-SLEX: black squares. 
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Figure 6: Error of numerical solutions at t = 12.5 for a circular orbit {R = 
1) as a function of scaled time step h (i.e., h = Xh). GR-IA: white discs 
(A = 0.8), GR-IA-LEX: gray discs (A = 1.1), GR-IA-SLEX: black discs 
(A = 1.5,1.8,1.5,1.7), GR-SYM: white squares (A = 1), GR-SYM-LEX: 
gray squares (A = 1.3,1.4,1.3,1.4), GR-SYM-SLEX: black squares (A = 
1.7,2.0,1.7,1.9). 
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